/* Bessel Y_1(x) function in double precision */

#include <math.h>

double dbesy1(double x)
{
    int k;
    double t, y, v, theta;
    static double a[16] = {
        -9.46499e-14, 2.141432796e-11, 
        -3.5977944347e-9, 4.3173541397185e-7, 
        -3.453883313621946e-5, 0.00165786399054057257, 
        -0.03978873577297383381, 0.31830988618379067154, 
        -5.571931463e-11, 8.93098810365e-9, 
        -9.9267351748541e-7, 7.164268731546281e-5, 
        -0.00295530533604595764, 0.05434868816050718549, 
        -0.1960570906462388432, -0.63661977236758134367
    };
    static double b[112] = {
        2.553685555e-11, 7.384080449e-11, 
        -2.9605446927e-10, -5.2081739106e-9, 
        -3.676554001376e-8, -1.532652533627e-7, 
        -1.3967947147386e-7, 3.88730181348638e-6, 
        4.052240780599446e-5, 2.5133005806772436e-4, 
        0.00109380629039004777, 0.00298157276348996763, 
        4.7188983238582062e-4, -0.19643830264444486066, 
        1.056916126e-11, 6.2160122556e-10, 
        3.85741197273e-9, 7.491469406e-9, 
        -6.306874425799e-8, -6.869833657799e-7, 
        -3.29785491054157e-6, -5.57893716372629e-6, 
        4.634569030169549e-5, 4.9086790937752985e-4, 
        0.00256488729380268318, 0.00822571374211646414, 
        0.01094512145046992425, -0.19159562910398869317, 
        -4.0211480479e-10, -2.29721816531e-9, 
        6.2931855684e-10, 6.509879184236e-8, 
        3.1589903600713e-7, -3.093448289001e-8, 
        -8.45532004628483e-6, -4.930719004112743e-5, 
        -9.865914153735679e-5, 4.7059803263198826e-4, 
        0.00472271367216891688, 0.01915518502418162917, 
        0.0372239136138415874, -0.1693323078896687811, 
        1.10728466147e-9, 2.05608441891e-9, 
        -2.536199750468e-8, -1.1751665325851e-7, 
        3.404534385077e-7, 4.36083390507752e-6, 
        8.34725302600388e-6, -7.154971448642437e-5, 
        -5.1956165332380432e-4, -0.00100811886858456098, 
        0.00436854257164588225, 0.03425644258611852849, 
        0.0907397117201493499, -0.10791597118914788657, 
        -1.17749768422e-9, 1.68403859797e-9, 
        3.471030171126e-8, -6.227816248757e-8, 
        -1.18405395385178e-6, 1.4478514453927e-7, 
        3.630744090017806e-5, 1.0481065165516454e-4, 
        -5.2483993430622472e-4, -0.00407014291708523893, 
        -0.00574311514714191566, 0.03534768318503458091, 
        0.16539360869993020163, 0.01986240425870529979, 
        4.008050552e-10, -3.36044652853e-9, 
        -9.2387130705e-10, 1.7840823064933e-7, 
        -2.2643988836534e-7, -8.08896250262948e-6, 
        -1.69403516892975e-6, 2.6617475073041589e-4, 
        7.2562705080786576e-4, -0.00399208199997748377, 
        -0.02399816920965700841, -0.0092557474455754244, 
        0.2008281402282066699, 0.21040564991987672576, 
        6.692975871e-11, 1.672202705e-10, 
        -1.433097981251e-8, 5.455077276306e-8, 
        1.04954046999028e-6, -3.39572219821681e-6, 
        -5.541577161092338e-5, 4.344113487321887e-5, 
        0.00184798594297387512, 0.0030104150846815045, 
        -0.02789739141907399514, -0.09421628519129373841, 
        0.09950212445529771576, 0.3749697583874640408, 
        4.528583204e-11, 6.6861092459e-10, 
        -9.2464816857e-9, -8.002058622489e-8, 
        8.7578676979688e-7, 6.28456150216647e-6, 
        -4.276239654834506e-5, -3.4630842915911896e-4, 
        8.9431202169021474e-4, 0.01086380107342193819, 
        0.00145563493724746238, -0.14193383817937981478, 
        -0.15148562296986494725, 0.35720232689066512346
    };
    static double c[126] = {
        3.201620498e-10, -1.77431944127e-9, 
        1.367520951871e-8, -6.171755212447e-8, 
        -4.3602779932091e-7, 1.79827173442894e-6, 
        4.990624965075363e-5, -1.8618453816525107e-4, 
        -0.00206578613620093339, 0.00640292258731756063, 
        0.04663627606907821596, -0.11400070635033154332, 
        -0.26159330264498333979, 0.30099732306965462346, 
        -5.9774708e-13, -3.8827648808e-10, 
        5.27269271362e-9, 2.971059814012e-8, 
        -5.2701406514027e-7, -3.20369808906187e-6, 
        4.475144442564492e-5, 1.6907703884518866e-4, 
        -0.00209937055264007564, -0.00491023437935628693, 
        0.0496717554070154514, 0.04195559522923170315, 
        -0.33516091307165838036, -0.02375823895638961835, 
        2.83380103e-12, -4.4715450235e-10, 
        1.769519961e-10, 6.05337984007e-8, 
        -2.888550814408e-8, -5.93780185020349e-6, 
        5.14027455732102e-6, 3.5667482807472357e-4, 
        -3.813272996260773e-4, -0.01158750857180657934, 
        0.01376562658346365929, 0.14388461069626297876, 
        -0.1310745466175553429, -0.27409127395927545296, 
        2.683697058e-11, -2.4302161327e-10, 
        -4.30151236236e-9, 3.557058903199e-8, 
        4.9328952592926e-7, -3.65732389065382e-6, 
        -3.642356546467712e-5, 2.3636175253326716e-4, 
        0.00154549935517002108, -0.00837084637397210913, 
        -0.02941202173717025197, 0.11713664042453548043, 
        0.15186375421302413107, -0.25912851048611625179, 
        2.914434604e-11, 1.5185121893e-10, 
        -4.88084673389e-9, -1.936460341189e-8, 
        5.7968735781303e-7, 1.58978381365925e-6, 
        -4.521858314881102e-5, -7.426575562367423e-5, 
        0.00206293925392279862, 0.00143934355204571641, 
        -0.04414786410004317489, -0.00317227451807568106, 
        0.273283773530321296, -0.02616867939853747003, 
        6.44550655e-12, 4.0249587279e-10, 
        -1.23789833341e-9, -5.58299757877e-8, 
        1.6975984093144e-7, 5.24003568623094e-6, 
        -1.539871212028104e-5, -3.0372172637669852e-4, 
        8.2284448843625174e-4, 0.00923614667011392611, 
        -0.02069442879928366384, -0.10834973445898588438, 
        0.14982326837249145873, 0.20317989938720766823, 
        -2.115691215e-11, 2.9752035314e-10, 
        3.34739862269e-9, -4.304625868398e-8, 
        -3.6733494871207e-7, 4.25328270524107e-6, 
        2.584872967637478e-5, -2.6244994622929809e-4, 
        -0.00102234403086242668, 0.00863185023199066935, 
        0.01816610124346617164, -0.11151661591132254182, 
        -0.0897879180557114995, 0.23370422835726857838, 
        -2.957977086e-11, -6.150420279e-11, 
        4.88089472193e-9, 6.24218182859e-9, 
        -5.6561995267846e-7, -3.2077350100271e-7, 
        4.27856984896483e-5, -5.3069698280036e-7, 
        -0.0018714153878560023, 7.3172821883575326e-4, 
        0.03832854206202077374, -0.01874044416229564036, 
        -0.2302705940514488064, 0.05794254714300082167, 
        -1.203179184e-11, -3.5509887475e-10, 
        2.13620090322e-9, 4.813181550803e-8, 
        -2.6815627351126e-7, -4.39191902440609e-6, 
        2.211921535574542e-5, 2.45935054805247e-4, 
        -0.00106158930833741898, -0.00722666167492552731, 
        0.02396424370522244155, 0.08278313570069734843, 
        -0.15890724632166919294, -0.15383825653750118007
    };
    static double d[52] = {
        -1.272346002224188092e-14, 3.370464692346669075e-13, 
        -1.144940314335484869e-11, 6.863141561083429745e-10, 
        -9.491933932960924159e-8, 5.301676561445687562e-5, 
        0.162867503967639974, -3.652982212914147794e-13, 
        1.151126750560028914e-11, -5.165585095674343486e-10, 
        4.657991250060549892e-8, -1.186794704692706504e-5, 
        0.01562499999999994026, 
        -8.713069680903981555e-15, 3.140780373478474935e-13, 
        -1.139089186076256597e-11, 6.862299023338785566e-10, 
        -9.491926788274594674e-8, 5.301676558106268323e-5, 
        0.162867503967646622, -2.792555727162752006e-13, 
        1.108650207651756807e-11, -5.156745588549830981e-10, 
        4.657894859077370979e-8, -1.186794650130550256e-5, 
        0.01562499999987299901, 
        -6.304859171204770696e-15, 2.857249044208791652e-13, 
        -1.124956921556753188e-11, 6.858482894906716661e-10, 
        -9.49186795351689846e-8, 5.301676509057781574e-5, 
        0.1628675039678191167, -2.185193490132496053e-13, 
        1.048820673697426074e-11, -5.132819367467680132e-10, 
        4.65740943737299422e-8, -1.186794150862988921e-5, 
        0.01562499999779270706, 
        -4.74041720979200985e-15, 2.578715253644144182e-13, 
        -1.104148898414138857e-11, 6.850134201626289183e-10, 
        -9.49167823417491964e-8, 5.301676277588728159e-5, 
        0.1628675039690033136, -1.75512205749384229e-13, 
        9.848723331445182397e-12, -5.094535425482245697e-10, 
        4.656255982268609304e-8, -1.186792402114394891e-5, 
        0.01562499998712198636
    };

    if (x < 0.85) {
        t = x * x;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * x;
        y = (((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15]) / x + 
            y * log(x);
    } else if (x < 4.15) {
        v = 4 / x;
        t = x - v;
        k = (int) (t + 4);
        t -= k - 3.5;
        k *= 14;
        y = (((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13]) * v;
    } else if (x < 12.5) {
        k = (int) x;
        t = x - (k + 0.5);
        k = 14 * (k - 4);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / x;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= -cos(x + theta);
    }
    return y;
}

